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Abstract 

Twisting and bending deformations are crucial to the biological functions 
of several microfilaments such as DNA molecules. Although continuum-rod 
models have emerged as efficient tools to describe the nonlinear dynamics 
of these deformations, a major roadblock in the continuum-mechanics-based 
description of microfilaments is the accurate modeling of the constitutive law, 
which follows from its atomistic-level structure and interactions. Since first- 
principle derivation of the constitutive law from atomistic-level structure and 
interactions is often impractical and so are direct experimental measurements 
due to the small length-scales, a natural alternative is to estimate the con- 
stitutive law from discrete-structure simulations such as molecular-dynamics 
(MD) simulations. In this paper, we describe a method based on a recently 
proposed two-step inverse approach [lj for estimating the constitutive law 
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using a static rod model and deformed configuration data generated from 
discrete-structure simulations. Furthermore, we illustrate the method on a 
filament with an artificial discrete-structure. We simulate its deformation 
in response to a prescribed loading using a multi-body dynamics (MBD) 
solver, Hyperworks MotionView. Using position data generated from the 
MBD solver, we first estimate the curvature of the filament, and subsequently 
use it in the described method to estimate the effective constitutive-law re- 
lationship between the restoring moment and curvature. Finally, we also 
illustrate how the estimated constitutive law can be validated under inde- 
pendent loading conditions. 

Keywords: Microstructure, constitutive behavior, rod theory, biological 
material, inverse methods. 



1. INTRODUCTION 

Microfilaments represent a wide range of natural or synthetic microscopic 
structures that are characterized by a large length-to-width aspect ratio and 
exhibit large bending and twisting deformations. Examples of microfila- 
ments include DNA molecules, microtubules, flagella, fimbriae and setae. 
The functions of these biological microfilaments are directly affected by their 
large nonlinear deformations and therefore it is important to understand and 
model the dynamics of these deformations. Although the ideas and methods 
discussed in this paper are applicable to almost all microfilaments, our mo- 
tivation stems from the fact that the biological functions of DNA molecules 
such as gene expression, transcription, and replication are significantly influ- 
enced by its long length-scale structural deformations such as looping 0, . 
Recent successes of continuum-rod models in describing biologically-relevant 
deformations of DNA molecules jl, [i| suggest that a continuum-rod model is 
an efficient tool to describe deformations of microfilaments. 

A continuum-rod model in general consists of dynamic equilibrium equa- 
tions and compatibility equations, which needs to be solved respecting a 
constitutive law. Although dynamic equilibrium equations and compatibil- 
ity equations, which we henceforth refer to as the rod model, remain the same 
for all microfilaments, the key distinguishing factor is the constitutive law. 
For each microfilament, the constitutive law, which describes the relationship 
between the deformation of the microfilament and the restoring forces and 
moments, follows from its atomistic structure and bond stiffnesses. 
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The rod model equations, that is, the equilibrium and compatibility equa- 
tions, reviewed in Section [21 can be recognized as four vector partial differ- 
ential equations in space s and time i |. In this paper, we assume the 
microfilament to be inextensible and unshearable. Under these assumptions, 
the four rod model equations involve five vector unknowns, namely, internal 
force f(s, t), internal moment q(s, t), curvature vector k(s, t), velocity v(s, t), 
and angular velocity u(s,t). The constitutive law serves as the fifth vector 
algebraic equation that captures the restoring effect of the internal force / 
and internal moment q on the curvature H. A general form of an elastic 
constitutive law applicable to an inextensible and unshearable rod model 
(discussed in Section [3]) is 

^ s) =0, (1) 

where tp is a vector function of its arguments and follow from atomistic-level 
interactions. While the rod model equations can be derived analytically 
from first principles, first-principle derivation of the constitutive law from 
atomistic-level interactions is often impractical. 

Alternatively, if all the arguments of ip, namely, internal force, inter- 
nal moment, and curvature vectors at each cross section, are experimentally 
measured, then well-established function-approximation tools can be used to 
estimate the constitutive law. Although such an approach will help circum- 
vent approximations due to modeling assumptions, it may introduce experi- 
mental artifacts. Nonetheless, while such measurements may be conceivable 
with acceptable accuracy for some laboratory-scale rods, the technology does 
not yet exist for sub-micron filaments. One of the many challenges at small 
length scales is that thermal fluctuations may dominate and thus corrupt any 
potential measurements. Some experiments exploit thermal fluctuations to 
measure persistence lengths and thus estimate average bending and torsional 
stiffnesses, which are dominant parameters in a linearized, homogeneous ap- 
proximation of equation (Tj[|); for example, refer to 0, S| for such experiments 
on DNA molecules. However, such simplistic constitutive laws are often in- 
adequate in modeling the deformation of microfilaments @, Eq|, EH, - 



A more practical alternative than using experimental measurements is 
to use data generated from discrete-structure simulations such as molecular 
dynamics (MD) simulations, which account for atomistic-level structure and 
interactions, and are accurate and computationally feasible for short length 



scales [13] . This alternative is already gaining attention for modeling con- 
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stitutive laws of several biological materials such as microtubules 14| and 



DNA [15j, although the attempts have been limited to estimation of only 
a few parameters, wherein the form of the constitutive law is presumed a 
priori, and is often presumed linear. Recently, a two-step inverse approach 
has been proposed |l| that can estimate directly the vector function ip from 
the discrete-structure simulations. This approach exploits state-space form 
of a static rod model, and can leverage several different scenarios of input 
and state measurements for the constitutive-law estimation. 

In this paper, we describe a method based on the two-step inverse ap- 
proach [H for a scenario, when the equilibrium positions of individual atoms 
computed from discrete-structure simulations under judiciously chosen load- 
ing conditions can provide the necessary information for estimating the ef- 
fective constitutive law. In this method, which is described in Section HJ we 
first estimate the curvature k from the position data of the deformed config- 
uration generated from discrete-structure simulations. This pre-processing; 



involves numerical differentiations and pose associated challenges (16| , which 



often necessitate smoothing of numerical data 0, Appendix B]. Next, follow- 
ing the step one of [if, we describe how the estimated curvature H can be 
used as a known quantity in the static rod model equations to solve for the 
internal restoring force / and internal restoring moment q in static equilib- 
rium of the deformed configuration. Finally, following the step two of [l[„ 
a least-squares fit can be used to estimate the functional relationship ip be- 
tween the estimated curvature k, the estimated force /, and the estimated 
moment q. 

In Section [5j we present results illustrating the applicability of the method 
on a filament with an artificial discrete-structure, which we simulate using 
a multi-body dynamics (MBD) solver. The artificial discrete-structure is so 
chosen that the expected constitutive law is homogeneous and has no depen- 
dence on internal force /. Furthermore, the artificial discrete-structure has a 
decoupled behavior in two-axes bending and torsion. So, it suffices to present 
the viability of the proposed method for planar bending about one principal 
bending axis. The simplicity of the artificial discrete-structure also allows 
us to offer a first-principle explanation of the estimated constitutive law as 
further corroboration. We, however, emphasize that the method described 
in Section H] is versatile, but illustrating its application to a variety of general 
cases is beyond the scope of this paper. 

To put our method in perspective, we begin by first reviewing a general- 
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purpose 3-D formulation of the rod model [6| in Section [2J which simulates 



the dynamics of large, nonlinear, deformations including intertwining [17 
In Section [31 we introduce a general, conceivable form of the constitutive law 
for an inextensible and unshearable microfilament. In Section HI we describe 
the method for the estimation of the constitutive law. The method exploits 
the static formulation of the rod model reviewed in Section [2J In Section El 
we present and discuss the results illustrating the estimation method on an 
artificial discrete-structure. We finally close by summarizing the conclusions 
in Section [6j 



2. The Rod Model 




Figure 1: Rod model of a microfilament in dynamic equiilibrium illustrating quantities of 
interest. 



In this section, we review the mathematical formulation of the rod model 
j6[ that involves dynamic equilibrium equations and compatibility conditions. 
In the rod model, the dynamics of the microfilament deformation follows 
from the rigid-body motion of its individual cross-sections. To track the 
rigid body motion of each cross-section, we fix at its mass center a reference 
frame cii(s,t), where the subscript i = 1,2,3, the independent variable s is 
the unstretched centerline coordinate and the independent variable t is time. 
As shown in Figure HJ we denote the position of the cross-section-fixed ref- 
erence frame by R(s,t). The gradient of R(s,t) along s, which we denote by 
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f(s, t), points along the centerline tangent t(s, t). Rod deformation in general 
involves two-axes-bending and torsion captured by a curvature vector H(s,t) 
as well as shear along the two flexure axes and extension (or compression) 
captured by f(s,t). In the stress-free (undeformed) state f(s,t) = t(s,t) and 
we denote the stress- free curvature by Ko(s). A change in the magnitude 
of r(s,t) captures extension (or compression), while change in its orienta- 
tion with respect to Oj(s, t) captures shear. The stress distribution across 
the cross-section results in a net internal (tensile and shear) force f(s,t) 
and (bending and torsional) moment q(s,t). The rigid body motion of the 
cross-section is described by its translational velocity v(s, t) and its angular 
velocity lo(s, t). 

The governing equations for the rod dynamics are 

df _ -* dv _ _ -> . . 

— + k x / = m— + mu x v — F, (2) 
as at 

— + nxq = l-— + uxl-uj + fxf-Q, (3) 
as at 

dv _ _ df . t . 

q^ + kxv = n+uxr, (4) 

where m(s) is the mass of the rod per unit length, the tensor I(s) is the 
mass moment of inertia per unit length, F(s,t) is the distributed external 
force per unit length, and Q(s,t) is the distributed external moment per 
unit length. The derivatives are all relative to the rotating, cross-section- 
fixed reference frame a,i(s,t). The dynamic equilibrium equations ()2]) and 
([3]) can be derived by applying Newton's second law to an infinitesimal rod 
element. The compatibility equations (jl]) and (jSJ) follow from the space- 
time continuity of the cross-section position R(s,t) and the cross-section 
orientation a,i(s,t) respectively 0, Appendix 2). A physical interpretation of 
these compatibility equations is that the gradient of the translational and 
angular velocities captures the time-rate of change of the rod's deformation. 
The rod model described above has been successfully exploited to of- 



fer promising insights on sequence-dependent DNA looping |4|, kinking [11 



intertwining 17| and buckling of microfibril arrays 18[|, albeit assuming sim 



plistic constitutive laws. Next, we discuss the constitutive law. 
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3. Constitutive Law for The Rod Model 



In this paper, we consider the rod model for slender microfilaments un- 
der low-tension applications that dominantly involve bending and twisting 
deformations. We assume that the microfilaments are inextensible and un- 
shearable. This assumption leads to the simplification f(s,t) = t(s,t), which 
is then a constant unit vector as seen from the cross-section-fixed reference 
frame aj(s,t) jl9[. Then, the rod model equations (EJ) - ([5]) constitute a set 
of four vector partial differential equations in five vector unknowns, namely 
translational velocity v(s,t), angular velocity u(s,t), internal force f(s,t), 
internal moment q(s,t), and curvature vector K,(s,t). The nonlinear rod 
model equations (0) - © can be solved numerically along with a fifth vec- 
tor equation, the constitutive law. A general form of an elastic constitutive 
law applicable to an inextensible and unshearable rod model is the algebraic 
constraint 

fU$J,s)=0. (6) 



The constitutive law (JH]) captures the restoring effects of internal force 
f(s,t) and internal moment q(s,t) on curvature vector K,(s,t) due to the 
stress-strain relationship integrated over the rod's cross section. It is worth- 
while to note that ignoring shear and extension does not ignore their kine- 
matic coupling with twisting and bending. This coupling is captured by 
internal force f(s,t) in the constitutive law (JH]). An example of a constitu- 
tive law with twist-extension coupling in inextensible rod model is presented 
in j^l, Section 6]. 

The elastic behavior of microfilaments and hence the form of their consti- 
tutive law follow from their atomistic- level interactions. In fact, there is no 
clear consensus on its functional form for DNA molecules and how it maps 
from the base-pair sequence. Although several experiments have provided 
estimates of average torsional and bending stiffness by measuring persistence 
lengths 0, 0], recent experimental observations [i[ and model- fitting 10, 12 



suggest that the constitutive law may be nonlinear, wherein the kink-ability 



of the molecule indicate non-convex stored energy functions [20]. 

Next, we present a general-purpose two-step computational method that 
can exploit the rod model equations along with data generated from discrete- 
structure simulations to estimate the vector function ijj. 
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4. Two-Step Estimation of Constitutive Law from Atomistic Struc- 
ture Simulations 

To estimate the effective constitutive law resulting from an atomistic pat- 
tern, we consider a cantilever microfilament of length L with several repeats 
of the atomistic pattern, and consider its deformation in static equilibrium. 
The cantilever is clamped at one end (at s = L) and body-fixed forces and 
moments are applied to the other end (at s — 0). This loading arrange- 
ment allows us to use a static formulation of the rod model, along with the 
data generated from discrete-structure simulations, in a two-step inverse ap- 
proach [l|, Table 1, Scenario 1] to estimate the effective constitutive law. For 
employing the two-step method, we write the static rod model equations in 
state-space form, which provides a general-purpose framework for a variety 
of estimation scenarios as elaborated in [H, Table 1]. 

4-1. Static Rod Model 

In static equilibrium, the rod model reduces to a time-independent for- 
mulation, that is, translational velocity v, angular velocity u and all partial 
derivatives with respect to t in equations (j2j) - (j3J) vanish. As a result, the 
compatibility equations (jlj) and (jSJ) are identically satisfied and the rod model 
reduces to two static equilibrium equations, namely 

~F, (7) 
fxf-Q (8) 

in three vector unknowns, namely internal force f(s), internal moment q(s) 
and curvature vector k(s). Since we assume that the rod is inextensible 
and unshearable, r = t, a known unit vector, which remains unchanged 
from reference (undeformed) configuration to deformed configuration as seen 
from the cross- sect ion-fixed reference frame di(s,t). We also recall that in 
equations ([7]) and (JS}, the derivatives are all relative to the cross-section- fixed 
reference frame di(s,t). 

4-2. State-Space Form of Static Rod Model 

To use the above static rod model equations, (J7J) and (JSJ) , in two-step 
estimation method, we refer them to cross- sect ion-fixed reference frame ct;(s) 



df 

ds 
dq 



+ KXf 

+ k x q 
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and write them in state-space form 



dx 

— = f{x,u) + w, se[0,L], (9) 

by defining the state x(s) £ R 6 as components of f(s) and q(s) along cij(s), 
input u(s) £ R 3 as components of k(s) along Oj(s) and disturbance iw(s) £ R 6 
as components of -F(s) and Q(s) along Oj(s). We recognize that a causal 
relationship need not exist between input u(s) and state x(s). In fact, the 
constitutive law (PQ) imposes an algebraic constraint between input u(s) and 
state x(s) 

$(x,u,s) = 0, $ £ R 6 x R 3 x R -> R 3 , (10) 

which for a repeating atomistic pattern will have a periodicity of A equal 
to the repeat interval and therefore $ (x, u,s) — $ (sc, it, s + A). Note that 
the state-space equation ([H]) is nonlinear, and moreover the state and input 
must satisfy the nonlinear algebraic constraint (fTUj) . We note that the state 
vector x(s) represents the force and the moment, the input u(s) represents 
the curvature vector and the disturbance w(s) represents distributed load on 
the cantilever. We also note that the state vector x(s) is prescribed at the 
free end of the cantilever (at s = 0), thus x(0) serves as initial conditions for 
the state-space equation (J9J) . 

If the constitutive law (JTUj) were known, the above state-space equations 
could be solved as a constrained initial value problem using a Differential 
Algebraic Equation (DAE) solver. Once the unknown curvature vector k(s) 
is solved, it can be further integrated to calculate the deformed shape of the 
rod respecting the position and orientation of the clamped end (at s = L) 



as summarized in Appendix A But since the constitutive law is unknown, 
we outline the following two-step method for estimating the constitutive law. 
We assume that the only measurements available from the discrete structure 
simulations are the reference (undeformed) and deformed configurations of 
the cantilever, distributed loading w(s) and the loading x(0) at the free end. 
Several other estimation scenarios based on measurement contingencies are 
discussed in U Table 1]. 

4-3. Two-Step Estimation Method 

In order to estimate the constitutive law ffTUj) . we first need to estimate its 
arguments, the state x(s) and the input u(s) from available measurements of 
a deformed state. Then, in the second step, we use function-approximation 
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Figure 2: Flowchart illustrating the components of the estimation of the constitutive law. 
The red arrows represent the path of information flow for a forward model simulation, while 
the blue arrows represent the path of information flow for inverse modeling to estimate 
the constitutive law. As seen from the chart, the inverse modeling framework uses data 
(or simulated data) and model information to estimate the constitutive law. For the 
results presented in this paper, the blocks shaded in green are implemented in Hypcrworks 
Motionvicw, while the blocks shaded in brown are implemented in MATLAB Simulink. 

tools such as tools based on least-square fitting to identify the functional 
relationship ( ITOj) . The first step, in general involves state-estimation and 
input-reconstruction methods depending on the available measurements [H, 
Table 1]. In this paper, we have position data as an available measurement. 
So, we derive the input u(s) (or equivalently, the curvature vector H(s)) 
from the position data using equation (JTBj) as explained next in Section |4.4p 
and then solve the state-space equation (Q for the state x(s). Thus in this 
scenario, we do not use state-estimation or input-reconstruction methods. 
Figure [2] illustrates the implementation of two-step estimation method in 
our scenario. 

We note that with the input u(s) known, equation (jHJ) becomes an un- 
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constrained initial value problem. We assume that accurate measurements of 
the initial condition x(0) (loading at the free end) and disturbance w(s) (dis- 
tributed loading) are available. Difficulties due to inaccurate/ incomplete 
measurements can be circumvented by leveraging partial measurements of 
state x(s) with unscented Kalman filter and/ or unscented unbiased minimum- 
variance (UUMV) filter as described in jlj]. 

4-4- Estimation of Curvature Vector from Position Data 

We use the position data of the reference (undeformed) and deformed 
configurations to estimate the input u(s) = {k(s)}. Here we introduced 
the notation {k(s)} to denote the 3x1 column matrix representing the 
components of the physical vector k(s) along a,i(s): 



f-1 




it-ax "J 








M 


K ■ a 2 ) 






, « • a 3 J 



(11) 



To derive {k(s)} from the position data, we first recall from p, Appendix 
2] that 



fl = -mn d2) 



where, the 3x3 square matrix 



e\ ■ a\ e 2 • a 1 e 3 • a\ 
ei ■ a 2 e 2 • a 2 e 3 ■ a 2 
h ■ a 3 e 2 ■ a 3 e 3 • a 3 



(13) 



represents the transformation from the inertial reference frame ij(s) (where 
j = 1, 2, 3) to the cross-section- fixed reference frame a)i(s) and 



r~i A 






~«3 


K 2 









-«a 








(14) 
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is 3x3 skew-symmetric matrix such that if vector c = k x b is the cross 
product of k with any vector b, then {c} = [k] j&j. The result ([T2"]) is 
derived in jsl, Appendix 2] from the Serret-Frenet formula 



and where the subscript ej denotes that the derivative is relative to the 
inertial reference frame ej. 

Next, let the square matrix [-R(s)] represent the rotation transforma- 
tiorfl of cross-section at s (or equivalently, of the cross-section-fixed frame 
Si(s)) from reference configuration to deformed configuration. Then [R] = 
[L re f] [Ldef] T and we derive [k(s)] from equation ffT2]) in terms of [R(s)] as 



The superscript T denotes matrix transpose, subscripts ref and def corre- 
spond to reference and deformed configurations respectively, and [k (s)] is 
skew-symmetric matrix corresponding to the stress-free (reference) curvature 
ko(s). Note that the stress- free curvature /?o(<s) depends on the choice of the 
cross-section-fixed reference frame cii(s) and satisfies equations ( |T2|) and ( TT5|) 
in reference configuration (with [L] = [L re f]). 

Finally, we recognize that if we estimate the rotation [-R(s)] from the ref- 
erence configuration to the deformed configuration, equation ffTB"]) provides 
the curvature vector k,(s) by calculating the skew-symmetric matrix [/c(s)], 
which has the three components of the curvature vector k(s) along a;(s) as 
per the equations fTTT]) and f TH|) . There are several conceivable ways of esti- 
mating the rotation [-R(s)] from the reference and the deformed configuration 
position-data. One approach is to select a point cloud (e.g. cluster of atoms) 
from the discrete structure around each cross-section and estimate the best 



1 Any change in the orientation of a cross-section from the undeformed configuration to 
the deformed configuration can be accomplished by a single Euler rotation about a unit 
vector u by angle 9. In terms of the single Euler rotation parameters, [R] = exp (9 [u]) = 
[I] + [u] sin 9+ [u] (I — cos 9), where [/] is the identity matrix and [u] is the skew-symmetric 
matrix corresponding to u. 




where recall that i = 1, 2, 3 



(15) 




(16) 
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fit of [-R(s)] for the point cloud from the reference configuration to the de- 
formed configuration. There are several commercially available tools (e.g. 
NX Imageware (UGS, Piano, TX)) to estimate the rigid body rotations of 
point clouds. 

An alternative approach to derive the curvature vector k(s) from position 
data is to use the Serret-Frenet formula (fT51) directly. This approach involves 
book-keeping of the cross-section- fixed reference frame fij(s) in deformed con- 
figuration from position data. We note that either approach involves numer- 
ical differentiations, which may potentially be ill-posed and amplify noise 



of the measured data [161 ] . The effects of noise in the two-step estimation 



method with appropriate unscented filters is analyzed in [1| 



5. Illustrative Results 

In this section, we illustrate the applicability of the two-step estimation 
method on a filament with an artificial discrete-structure with planar bend- 
ing. We first describe the artificial discrete-structure in Sub-section 15.11 
Next, in Sub-section 15.2} we exploit some simplicities in the expected form of 
the constitutive law that are physically obvious from the discrete structure. 
In Sub-section 15.31 we describe the state-space form of the reduced-order 
static rod model in the plane of bending. Finally, we present and discuss the 
estimation and validation results in Sub-sections 15.41 and 15.51 respectively. 

5.1. An artificial discrete- structure 

We construct and simulate the discrete structure in a commercial, multi- 
body-dynamics software, Hyperworks Motion View. We follow the default 
system of units in Hyperworks, which is [Newton, Millimeter, Kilogram, Sec- 
ond] . In the discrete structure, we have thirty repeats of a cubic arrangement 
of unit point masses connected by linear springs of unit stiffness. Figure [3] 
shows two orthogonal views of a representative portion (five repeats) of the 
discrete structure in undeformed (reference) state. The springs are along 
each edge of the cube as well as along diagonals of two opposite faces as 
shown. Each edge of the cube is of unit length, so, the length of the entire 
filament L = 30 mm. 

5.2. Simplified constitutive law with in-plane bending 

Before estimating the effective constitutive law, we first exploit the sim- 
plicities in its form that are physically obvious from the discrete structure. 
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Figure 3: Discrete structure of the filament. Two orthogonal views are shown, one along 
z-axis (top) and the other into the z-axis (bottom). Unit point masses are arranged in 
a cubic arrangement and five repeats of the cubic arrangement are shown. Each edge of 
the cubic cell has unit length and has linear springs of unit stiffness connecting the point 
masses in the undeformed (reference) configuration. In addition, the diagonal springs are 
placed on two opposite faces of the cube shown in the view along z-axis (top), but not on 
the other four faces of the cube. 

Since the discrete structure is uniform, we expect the effective constitutive 
law also to be uniform (homogeneous) along the length. In other words, there 
is no explicit dependence on s in equation (JT]). In addition, the constitutive- 
law equation ([ID also has no explicit dependence on the internal force f(s) 
for the chosen discrete structure 0. So, the equation flT| that needs to be 
estimated, simplifies to the following form: 

f{*,<f) = 0- (17) 

Finally, we also note that due to the symmetries of the discrete-structure, 
the constitutive law is decoupled in the principal directions of bending and 
torsion. Choosing the cross-section-fixed reference frame a,i(s) along the prin- 
cipal directions, the vector equation (|T7|) simplifies to three scalar equations 

ipi(Ki,qi) = 0, (18) 

where the subscript i = 1,2, or 3 corresponds to each of the principal axes 
along ai(s). Specifically, we let fii(s), a 2 (s), and a 3 (s) correspond to the 
first bending, second bending, and torsion axes, respectively, as illustrated 
along with a result in Figure HI The unit vector a 3 (s) points along the 



2 Although this simplicity may not be readily obvious, all the simplicities that we impose 
in this section will become plausible anyway with validation of the estimated constitutive 
law in Section [5.51 
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outward normal to the cross-section. The unit vector ai(s) is parallel to the 
faces of the cube containing diagonal springs, while a 2 (s) is parallel to the 
faces not containing diagonal springs To present the viability of the two- 
step estimation method, we illustrate it for planar bending about only one 
principal bending axis, which we choose to be along a2(s)-axis. 

5. 3. In-plane static rod model in state-space form 

In-plane bending about the ct2-axis occurs when a shear force is applied 
in the ai-direction and/ or bending moment is applied about the a 2 -axis. 
Referring all the physical vectors to the cross-section-fixed reference frame 
ctj(s), we note that the first and third components of k(s) and q(s) and the 
second component of f(s) vanish for the planar bending. Furthermore, we 
will have no distributed load. In this case the governing vector equations ([7]) 
and ([H]) of the static rod model further reduce to the following three scalar 
nonlinear differential equations 

dq 2 
ds 
dfi 
ds 
dfs 
ds 

with the constitutive law to be estimated 

V> 2 (K2,<?2) = 0. (22) 



-fl, 


(19) 


-/3K2, 


(20) 


/l«2, 


(21) 



Equations ( fT9TT2~TT) are in the state-space form ([9]), as described in [2l|], with 











X = < 




H 













and u = «2- 



(23) 



5.4- Estimation of the constitutive law 

For estimation of the effective constitutive law, we deform the cantilevered 
discrete-structure in Hyperworks MotionSolve by exerting a shear force /i(0) = 



3 Note that with this choice of reference frame, the stress-free curvature vector kq(s) = 
since the discrete structure is straight in undeformed (reference) state. This fact can be 
verified by satisfying Serret-Frenet formula (fT"5j) in the reference state. 
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Figure 4: Deformed shape of the discrete-structure cantilever in static equilibrium. The 
free end has a prescribed shear force, but has no tension and no bending moment. Since 
there arc thirty cubic cells of 1 mm edge each, L = 30 mm. The cross-section-fixid reference 
frame di(s) is also shown. The unit vector 03 (s) points along the outward normal to the 
cross-section. The unit vector ai(s) is chosen in the plane of bending, while 0,2 (s) is parallel 
to the axis of bending. 



—0.004 N at the free end (at s = 0) as shown in Figure HI We simulate the 
static equilibrium configuration of the discrete structure using the Force Im- 
balance Method- Type D for the static analysis in Hyperworks MotionSolve. 
From the deformed configuration, we estimate the curvature K2{s) as out- 
lined in Section H31 Next, to estimate the internal moment q2 (s), we use the 
estimated curvature k 2 (s) as a known input u(s) in the state-space equations 
( fl9|) - (121~|) of the in-plane static rod model. We solve the state-space equa- 
tions ( II 9p - (12 ip for x(s), which contains 92 (s), in MATLAB Simulink with 
the initial conditions: 
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(rad/mm) 

Figure 5: Estimated constitutive law. 



f /l(0) 1 A f -°- 004 N ] 

x(0) = /a(0) I = ON I (24) 

[ g 2 (0) J [ N-mm J 

Finally, to estimate the functional relationship (122]) between and #2 
we express the function, ip2, as a sinusoidal-basis-function expansion, where 
the unknown coefficients of this expansion are then found by standard least- 
squares-fitting. Figure [5] shows the estimated constitutive relationship be- 
tween the bending moment and curvature. 

5. 5. Validation of the constitutive law 

To validate the estimated constitutive-law, one may test it in several 
different loading environments (static or dynamic) and evaluate if it predicts 
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(mm) 



Figure 6: Deformed shape of the discrete-structure cantilever in static equilibrium (left). 
The free end has prescribed shear force, tension as well as bending moment. Center- 
line predicted from discrete-structure simulation in Hypcrworks (dots) and continuum-rod 
simulation (solid curve) employing the estimated constitutive-law (right). 

the same deformation behavior from the rod model as does the discrete- 
structure simulation. Here we present one spot test in static equilibrium. The 
loading conditions for the spot test are shown in Figure El on the deformed 
equilibrium of the discrete structure (on the left). In this case we apply shear 
force, compressive force as well as bending moment at the free end. We apply 
these loading conditions in the Simulink state-space-model of the continuum 
rod employing the estimated constitutive- law shown in Figure El We set it up 
as DAE in Simulink with constitutive law as the known, algebraic constraint. 
Figure M shows (on the right) the comparison of the centerline predicted from 
the discrete-structure simulation in Hyperworks (dots) and the continuum- 
rod simulation (solid curve) employing the estimated constitutive-law. The 
closely matching shape validates the estimated constitutive- law. 
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5.6. Analysis and discussion of results 

The close agreement obtained in Section |5"31 of the centerline deformations 
of the discrete-lattice structure and the continuum rod-model in the spot-test 
loading not only illustrates the viability of the two-step estimation-method, 
but also corroborates several assumptions that we made along the way. In 
particular, we assumed inextensibility and unshearability in the rod model 
equations ( 1T9"|) - (12Tp . while the discrete-structure simulations are inherently 
free from such assumptions. The slenderness ratio (length/thickness) of the 
discrete structure is 30. Therefore the results corroborate the assumptions 
of inextensibility and unshearability in Kirchhoff rod theory for slender fil- 
aments. Furthermore, the close agreement corroborates the simplifications 
that we imposed a priori on the constitutive law in Section 15.21 The agree- 
ment is also noteworthy when considering that each simulation is run in- 
dependently with different softwares and solvers, and therefore serve as a 
validation study for the softwares and solvers as well. The accuracy of the 
results obtained using commercial softwares and standard numerical meth- 
ods is promising for the ease of implementing the method to a variety of 
applications involving filament-type deformable structures. 

We recall here that the first-principle derivation of the constitutive law 
from atomistic-level interactions is often impractical in general. However, for 
the discrete structure that we chose, we indeed can corroborate the estimated 
constitutive law from an approximate first-principle derivation. Consider the 
spring forces in the deformed discrete structure in Figure [7] showing one 
cubic cell. If R is the radius of curvature of the neutral surface and 9 is the 
bend angle, R9 = 1 mm and k 2 — 1/R = 6 /mm. The restoring moment 
q 2 is related to the restoring forces developed in the springs on the top and 
bottom faces of the cube. The top two springs in the plane of bending are 
stretched, while the bottom two springs are compressed. The restoring force 
/ developed in each of the four springs is 9/2 N and therefore, the resultant 
moment q 2 = 9 N-mm. Eliminating 9 between q 2 and k 2 , we get q 2 = (1 
N-mm 2 ) k 2 , which corroborates the nearly linear relationship between q 2 and 
k 2 observed in Figure EJ 

6. Conclusions 

The continuum-rod model has emerged as an efficient tool to describe the 
long length and time scale dynamics of nonlinear deformations of microfila- 
ments and sem-flexible biofilaments such as DNA. However, the usefulness 



19 




Figure 7: A schematic to derive the constitutive law ifa («2, (I2) = from the restoring 
spring-forces developed in the discrete structure. The neutral surface is approximately 
half way between the top and bottom surfaces. 

of the rod model is severely limited by the lack of knowledge of an accurate 
constitutive law. Discrete-structure simulations, such as MD simulations of 
biomolecules such as DNA for small length and time scales can provide ample 
data for the accurate estimation of an effective constitutive law. 

In this paper, we described a method based on a recently proposed two- 
step inverse approach [l[ to estimate the constitutive law using discrete- 
structure simulations. We illustrated the applicability of the method on a 
simple representative example. In particular, we considered a cantilever fil- 
ament with an artificial discrete-structure. We simulated its deformation in 
response to a prescribed loading using a multi-body dynamics (MBD) solver. 
We then estimated an effective constitutive law from the deformation data 
using the two-step method. Finally, we illustrated how the estimated consti- 
tutive law can be validated by employing it in the continuum-rod model and 
comparing the simulation results with those of discrete-structure simulations 
under an independent set of cantilever loading conditions. The method pre- 
sented in this paper can leverage MD simulations of the discrete atomistic 
structure of microfilaments to estimate their constitutive laws. 
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Appendix A. Shape from Curvature 

Appendix A.l. Three-dimensional deformation 

Here, we summarize how to solve for the centerline shape R(s) from the 
curvature components {k(s)} for an inextensible and unshearable rod. Recall 
from Section [2] that 

= r, (A.l) 

where the subscript ej denotes that the derivative is relative to the inertial 
reference frame ej. For an inextensible and unshearable rod, the position 
gradient f(s) = t(s). So, we can integrate equation ( lA.ip to get R(s) if we 
find 





(A.2) 

the components of the unit tangent vector t(s) along ij. But 



{% = [Lf{i}, (A.3) 

where [L(s)\ is the transformation matrix defined in equation flT3|) . the su- 
perscript T denotes the transpose and the 3x1 matrix 
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{i- di "I 
t-a 2 \ (A.4) 
i-a 3 J 

represents the components of i(s) along hi. We know {£} a priori from the 
choice of hi. The transformation matrix [L(s)} can be solved from {k,(s)} by 
integrating equation (fl2l) . We use the following approximation to integrate 
equation ( fl~2l : 



[L(s + <5s)] « exp 



(A.5) 



where 



is the skew-symmetric matrix corresponding to 



A 



s+Ss 

{ft} ds (A.6) 



such that if vector c = 9 x is the cross product of # with any vector b, then 



{c} 







We evaluate the integral in equation ( 1A.6I) using trapezoidal rule, inte 



grate equation ( 1X21) to solve for [L(s)] using numerical approximation ( 1A.5I) . 
solve for {t(s)}„ from equation f ]A.3j) and finally recognizing r(s) = t(s), 

integrate equation flA.lj) using trapezoidal rule to get the centerline curve 
R{s) of the deformed rod. 

Appendix A. 2. In-plane bending 

For in-plane bending presented in Section [51 the shape integration from 
the curvature k 2 is simpler than the above algorithm for the general 3-D case. 
If the plane of bending is designated as x — y plane and if 0(s) is the angle 
that centerline tangent makes with the x-axis, then 




(A.7) 



Furthermore, if the x and y components of R{s) are denoted by R x {s) and 
R y (s) respectively, then 
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= cos0 and ^^- = sin0. (A.i 
ds ds 



Thus, we can integrate equation (1A.7j) using trapezoidal rule to get 0(s 



and then integrate equation (1A.8I) to get R x (s) and R y (s) respecting the 
boundary conditions at the clamped end (at s — L) of the cantilever. In 
particular, we can integrate along the increasing s-direction starting from the 
free end (s = 0) with 0(0), -R^O) and R y (0) (the constants of integration) as 
unknowns and find their values satisfying the clamped values of 4>(L), R X (L) 
and R y (L). 
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